



load(paste0(datDir,"estimDat_minPpl_all_5.RData"))

estimDatUse = estimDat[mult_admit == 3 & age06 <= ageCut06 ,]
## Developed in fl_compute_deltas.R
load(paste0(datDir,"deltaDat_minPpl_all_5.RData"))

source('./SwitchCosts/stata_results_switch.R')
DISTANCE = coefDistList$main
SWITCHING.COST =  coefSwitchList$main


data = grpShrDatXi[is.finite(xiHat),][.id == "main",]
setnames(data,"nGrp","nZip")

dataMkt = data

dataMkt[,xiHat := deltaHat - DISTANCE * (time_current)]
dataMkt  = dataMkt[,list(hosp_id_alt,pat_zip,xiHat,time_current,deltaHat)]
## Add Outside Option
dataMkt  = rbind(dataMkt,as.data.table(list(hosp_id_alt = "0",xiHat = 0, time_current = 0, deltaHat = 0, pat_zip = dataMkt[,unique(pat_zip)])))

dataMktMerge = dataMkt[,list(pat_zip,hosp_id_alt,deltaHat,xiHat)]
dataMktMerge[,hosp_id_alt := as.numeric(hosp_id_alt)]
dataMerged = copy(merge(estimDatUse,dataMktMerge, by = c("pat_zip","hosp_id_alt")))


setorder(dataMerged, maskssn, hosp_id_alt,birthNum)




doSimulation = function(ITER,data,YEAR.START,YEAR.END,DISTANCE, SWITCHING.COST,sdXi, sdPsi, correlPsi) {
 
  dataSimul = copy(data)
  dataSimul[birthNum == 1,xiHatNew := xiHat + rnorm(n = length(dataSimul[birthNum == 1,xiHat]),mean = 0 ,sd = sdXi)]
  dataSimul[,xiHatNew := mean(xiHatNew, na.rm = TRUE), by = c("maskssn","hosp_id_alt")]
  dataSimul[birthNum == 1,psiHat := rnorm(n = length(dataSimul[birthNum == 1,xiHat]),mean = 0 ,sd = sdPsi/sqrt(1 - correlPsi^2))]
  dataSimul[, psiHatOld := c(NA,head(psiHat,-1)), by = list(maskssn,hosp_id_alt)]
  dataSimul[birthNum == 2, psiHat := psiHatOld * correlPsi +  rnorm(n = length(dataSimul[birthNum == 2,xiHat]),mean = 0 ,sd = sdPsi)]
  dataSimul[, psiHatOld := c(NA,head(psiHat,-1)), by = list(maskssn,hosp_id_alt)]
  dataSimul[birthNum == 3, psiHat := psiHatOld * correlPsi +  rnorm(n = length(dataSimul[birthNum == 3,xiHat]),mean = 0 ,sd = sdPsi)]
  dataSimul[,rand := runif(1), by = c("maskssn","birthNum")]
  for(YEAR in YEAR.START:YEAR.END) {
    dataSimul[yearDat == YEAR, detUtil := xiHatNew + psiHat + DISTANCE * (time_current) + SWITCHING.COST * (chosenPrev == TRUE)]
    ## Outside option has no utility, no switching benefit
    dataSimul[yearDat == YEAR & hosp_id_alt == 0, detUtil := 0]
    
    dataSimul[yearDat == YEAR, hospShare := exp(detUtil)/sum(exp(detUtil)), by = c("maskssn","birthNum")]
    dataSimul[yearDat == YEAR, hospShareCumSum := cumsum(hospShare), by = c("maskssn","birthNum")]
    dataSimul[yearDat == YEAR, obs_choice := (rand<= hospShareCumSum & rand > (hospShareCumSum - hospShare))]
    setkey(dataSimul,maskssn,hosp_id_alt,yearDat)
    ## Overwrites some of the old chosenPrev
    dataSimul[, chosenPrevNew:= c(NA,head(obs_choice,-1)), by = list(maskssn,hosp_id_alt)]
    dataSimul[!(is.na(chosenPrevNew) | hosp_id_alt == "0"), chosenPrev := chosenPrevNew]
  }
  
  dataSimul[mult_admit > 1 & birthNum <= 2 & hosp_id_alt != 0, nHospChosen := sum(obs_choice), by=list(maskssn, hosp_id_alt)]  ##Number of hospitals chosen in first two births
  dataSimul[mult_admit > 1,hospSwitchFirst2 := (sum(nHospChosen == 1, na.rm=TRUE) == 4), by = 'maskssn'] ##People who went to different hospitals for first two births, where both were in choice set both times
  setkey(dataSimul, maskssn, hosp_id_alt)

  dataSimul[mult_admit > 1,chosenFirstHospTemp := (obs_choice == TRUE & birthNum == 1)]
 
  dataSimul[mult_admit > 1,chosenSecHospTemp := (obs_choice==TRUE & birthNum == 2)]
  

  varNames <- c('chosenFirstHosp','chosenSecHosp')
  varNamesTemp <- paste0(varNames, 'Temp')
  for (nm in varNames){
    dataSimul[mult_admit > 1, (nm) := as.integer(max(get(paste0(nm,'Temp')), na.rm=TRUE)),by = list(maskssn,hosp_id_alt)]
  }
  dataSimul[birthNum == 3 & mult_admit > 1,firstHosp3Temp := as.integer(chosenFirstHosp == 1 & obs_choice == TRUE)]

  dataSimul[birthNum == 3 & mult_admit > 1,secondHosp3Temp := as.integer(chosenSecHosp == 1 & obs_choice == TRUE)]

  dataSimul[birthNum == 2 & mult_admit > 1,noMove2Temp := as.integer(pat_zip == pat_zip_prev)]

  dataSimul[birthNum == 3 & mult_admit > 1,noMove3Temp := as.integer(pat_zip==pat_zip_prev)]
  
  
  
  varNames <- c('firstHosp3','secondHosp3','noMove2','noMove3')
  varNamesTemp <- paste0(varNames, 'Temp')
  for (nm in varNames){

    dataSimul[mult_admit > 1, (nm) := as.integer(max(get(paste0(nm,'Temp')), na.rm=TRUE)),by = list(maskssn)]
  }
  

  for (nm in grep('Temp',names(dataSimul), value = TRUE))  dataSimul[,(nm) := NULL] 
  
  

  chambDatMelt <- melt(dataSimul[hospSwitchFirst2 == TRUE & 
                                       (chosenFirstHosp == 1 | chosenSecHosp == 1) & birthNum <=3,
                                     list(maskssn,birthNum,age06, time_current, chosenFirstHosp,noMove3,
                                          firstHosp3,secondHosp3,
                                          pat_zip,hosp_id_alt)],
                       id.vars=c('maskssn','birthNum', 'age06', 
                                 'chosenFirstHosp','noMove3','firstHosp3','secondHosp3'))
  
  chambDat <- data.table(dcast(chambDatMelt,
                               age06 + maskssn+noMove3+  firstHosp3+secondHosp3  ~
                                 variable+birthNum+chosenFirstHosp))
  
  chambDat[,one := 1]
  chambDat[,time_current_diff:=(time_current_2_0-time_current_2_1)-(time_current_1_0-time_current_1_1)]
  
  honoreDatUse <- subset(chambDat, age06 <= ageCut06 & noMove3==1 & hosp_id_alt_1_0 != hosp_id_alt_1_1 )
  honoreDatUse[,secondHospMinFirstHosp3 := (secondHosp3==1)-(firstHosp3==1)]
  
  honoreEst = coef(glm(one~ secondHospMinFirstHosp3 + time_current_diff- 1, family = 'binomial', data = honoreDatUse))
  return(honoreEst)
}

doSimulationWrap = function(ITER,data,DISTANCE,SWITCHING.COST,YEAR.START,YEAR.END,sdXi, sdPsi, correlPsi) {
  dataSimulList= list()

    set.seed(123457 + ITER)

    results = doSimulation(ITER = ITER,data = data, DISTANCE = DISTANCE,SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END,sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi)
  
return(results)

}


YEAR.START = 2006
YEAR.END = 2014


sdXi = 0
sdPsi = 0
correlPsi = 0

time.start = proc.time()
ITER.LIMIT = 100
simulTest = NULL
for(ITER in 1:ITER.LIMIT) {
  simulTest =  rbind(simulTest,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start


saveRDS(simulTest,paste0(resultsDir,"/simulMC_Test.rds"))


sdXi = 0.5
sdPsi = 0
correlPsi = 0

time.start = proc.time()
ITER.LIMIT = 100
simulTest2 = NULL
for(ITER in 1:ITER.LIMIT) {
  simulTest2 =  rbind(simulTest2,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start


saveRDS(simulTest2,paste0(resultsDir,"/simulMC_Test2.rds"))



sdXi = 0.5
correlPsi = 0.1
sdPsi = 0.1*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul = NULL
  for(ITER in 1:ITER.LIMIT) {
simul =  rbind(simul,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start


saveRDS(simul,paste0(resultsDir,"/simulMC_1.rds"))



sdXi = 0.5
correlPsi = 0.1
sdPsi = 0.5*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul2 = NULL
for(ITER in 1:ITER.LIMIT) {
  simul2 =  rbind(simul2,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start


saveRDS(simul2,paste0(resultsDir,"/simulMC_2.rds"))


sdXi = 0.5
correlPsi = 0.5
sdPsi = 0.5*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul3 = NULL
for(ITER in 1:ITER.LIMIT) {
  simul3 =  rbind(simul3,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start

saveRDS(simul3,paste0(resultsDir,"/simulMC_3.rds"))

sdXi = 0.5
correlPsi = 0.5
sdPsi = 0.1*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul4 = NULL
for(ITER in 1:ITER.LIMIT) {
  simul4 =  rbind(simul4,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start

saveRDS(simul4,paste0(resultsDir,"/simulMC_4.rds"))

sdXi = 0.5
correlPsi = 0.5
sdPsi = 1*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul5 = NULL
for(ITER in 1:ITER.LIMIT) {
  simul5 =  rbind(simul5,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start

saveRDS(simul5,paste0(resultsDir,"/simulMC_5.rds"))

sdXi = 0.5
correlPsi = 0.9
sdPsi = 1*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul6 = NULL
for(ITER in 1:ITER.LIMIT) {
  simul6 =  rbind(simul6,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start

saveRDS(simul6,paste0(resultsDir,"/simulMC_6.rds"))

sdXi = 0.5
correlPsi = 0.1
sdPsi = 1*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul7 = NULL
for(ITER in 1:ITER.LIMIT) {
  simul7 =  rbind(simul7,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start


saveRDS(simul7,paste0(resultsDir,"/simulMC_7.rds"))


sdXi = 0.5
correlPsi = 0.9
sdPsi = 0.5*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul8 = NULL
for(ITER in 1:ITER.LIMIT) {
  simul8 =  rbind(simul8,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start


saveRDS(simul8,paste0(resultsDir,"/simulMC_8.rds"))


sdXi = 0.5
correlPsi = 0.9
sdPsi = 0.1*sqrt(1 - correlPsi^2)


time.start = proc.time()
ITER.LIMIT = 100
simul9 = NULL
for(ITER in 1:ITER.LIMIT) {
  simul9 =  rbind(simul9,doSimulationWrap(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END, sdXi = sdXi, sdPsi = sdPsi, correlPsi = correlPsi))
}
time.finish = proc.time()
time.finish - time.start


saveRDS(simul9,paste0(resultsDir,"/simulMC_9.rds"))

simulResults = as.data.table(rbind(cbind(as.data.table(simulTest[1:ITER.LIMIT,]),data.table(sdXi = 0,sdPsi = 0, correlPsi = 0)),
                                   cbind(as.data.table(simulTest2[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 0, correlPsi = 0)),
                                   cbind(as.data.table(simul[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 0.1, correlPsi = 0.1)),
                                   cbind(as.data.table(simul2[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 0.5, correlPsi = 0.1)),
                                   cbind(as.data.table(simul3[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 0.5, correlPsi = 0.5)),
                                   cbind(as.data.table(simul4[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 0.1, correlPsi = 0.5)),
                                   cbind(as.data.table(simul5[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 1, correlPsi = 0.5)),
                                   cbind(as.data.table(simul6[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 1, correlPsi = 0.9)),
                                   cbind(as.data.table(simul7[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 1, correlPsi = 0.1)),
                                   cbind(as.data.table(simul8[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 0.5, correlPsi = 0.9)),
                                   cbind(as.data.table(simul9[1:ITER.LIMIT,]),data.table(sdXi = 0.5,sdPsi = 0.1, correlPsi = 0.9))))
setnames(simulResults,c("secondHospMinFirstHosp3","time_current_diff"),c("sc","time"))                                  
                                   
simulResults = data.table(simulResults)

ggplot(data = simulResults[correlPsi > 0,],aes(x = as.factor(sdPsi),y = sc)) + geom_boxplot() +
  geom_hline(yintercept = SWITCHING.COST, color = "red") + theme_bw() + facet_grid( .~ correlPsi) +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  labs(x = "SD of Psi", y = "Switching Cost") + scale_y_continuous(breaks = c(0.7, 0.8, 0.9, 1,1.1)) + expand_limits(y = c(0.65, 1.15))
ggsave(filename = paste0(resultsDir,"/Switch/MC_Psi_SC","_Pres",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/MC_Psi_SC","_Paper",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")

ggplot(data = simulResults[correlPsi > 0,],aes(x = as.factor(sdPsi),y = -1*time)) + geom_boxplot() +
  geom_hline(yintercept = -1*DISTANCE, color = "red") + theme_bw() + facet_grid( .~ correlPsi) +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  labs(x = "SD of Psi", y = "Distance") + scale_y_continuous(breaks = 1*c(0.03,0.04, 0.05, 0.06, 0.07,0.08,0.09,0.10,0.11)) + expand_limits(y = 1* c(0.03, 0.11))

ggsave(filename = paste0(resultsDir,"/Switch/MC_Psi_Distance","_Pres",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/MC_Psi_Distance","_Paper",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")

ggplot(data = simulResults[correlPsi == 0,],aes(x = as.factor(sdXi),y = sc)) + geom_boxplot() +
  geom_hline(yintercept = SWITCHING.COST, color = "red") + theme_bw() +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  labs(x = "SD of Xi", y = "Switching Cost") + scale_y_continuous(breaks = c(0.7, 0.8, 0.9, 1,1.1)) + expand_limits(y = c(0.65, 1.15))

ggsave(filename = paste0(resultsDir,"/Switch/MC_Xi_SC","_Pres",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/MC_Xi_SC","_Paper",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")

ggplot(data = simulResults[correlPsi == 0,],aes(x = as.factor(sdXi),y = -1*time)) + geom_boxplot() +
  geom_hline(yintercept = -1*DISTANCE, color = "red") + theme_bw() +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  labs(x = "SD of Xi", y = "Distance") + scale_y_continuous(breaks = 1*c(0.03,0.04, 0.05, 0.06, 0.07,0.08,0.09,0.10,0.11)) + expand_limits(y = 1* c(0.03, 0.11))


ggsave(filename = paste0(resultsDir,"/Switch/MC_Xi_Distance","_Pres",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/MC_Xi_Distance","_Paper",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")
